function generateMfuncResidualEquation_file(allModelParams,x,y,xp,yp,nx1,ne,h11,strucOfArrays,positionFile)
% This function generates an M function which is the objective function
% when solving the model by projection
% The results are saved in the file ResidualEquationProj.m

% We start deleting the old version of the file - if it exists
nameOfFunction = 'ResidualEquationProj.m';
saveFile = [positionFile,'\',nameOfFunction];
if exist(saveFile,'file') > 0
    delete(saveFile)
end

text = ['function [res,resEuler,dF] = ',nameOfFunction(1:end-2),'(thetaVec,setup)'];
dlmwrite(saveFile,text,'-append','delimiter','');
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '%% Getting some information from setup ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'scaleX_in_g = setup.scaleX_in_g;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'x_cu        = setup.xGrid;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'numX        = size(x_cu,1);';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'nNodes      = setup.n_nodes;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'wNodes      = transpose(setup.weight_nodes);';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'epsNodes    = setup.epsi_nodes;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'params      = setup.params;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'appOrder    = setup.appOrder;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'nx          = setup.nx;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'ny          = setup.ny;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'nx1         = setup.nx1;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'gSS         = setup.gSS;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'hSS         = setup.hSS;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'thetaProj   = reshape(thetaVec,setup.nDim1Theta,setup.nDim2Theta);';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'g0          = thetaProj(:,1); ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'gx          = thetaProj(:,2:1+nx); ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'if appOrder > 1 ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '    gxxV    = thetaProj(:,nx+2:1+nx+nx*(nx+1)/2); ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'end ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'if appOrder > 2 ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '    gxxxV   = thetaProj(:,1+nx+nx*(nx+1)/2+1:1+nx+nx*(nx+1)/2+nx*(nx+1)*(nx+2)/6); ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'end ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '%Unfold params ';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:length(allModelParams)
    text = [allModelParams{i}, '= params.',allModelParams{i},';'];
    dlmwrite(saveFile,text,'-append','delimiter','');
end 
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '%% Building the residuals ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '% Getting the states at time t ';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:length(x)
     text = [char(x(i)),' = hSS(',num2str(i),',1) + repmat(x_cu(:,', num2str(i),'),1,nNodes);'];
     dlmwrite(saveFile,text,'-append','delimiter','');
end
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '% Getting the controls at time t and part of the states at time t+1  ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '% yVec_cu: numX * ny ';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:length(x)
    text = ['x',num2str(i),'_cu ',' = x_cu(:,',num2str(i),')/scaleX_in_g(',num2str(i),',1);'];
    dlmwrite(saveFile,text,'-append','delimiter','');
end
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:length(y)
    name = [char(y(1,i)),' = ','g0(',num2str(i),',1)'];
    for j=1:length(x)
        name = [name,['+gx(',num2str(i),',',num2str(j),').*x',num2str(j),'_cu']];
    end
    text = [name,';'];
    dlmwrite(saveFile,text,'-append','delimiter','');
end
text = 'if appOrder > 1';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:length(y)
    name = [char(y(1,i)),' = ',char(y(1,i))];
    idx = 0;
    for alfa1=1:length(x)
        for alfa2=alfa1:length(x)
            idx = idx + 1;
            if alfa1 == alfa2
                name = [name,['+gxxV(',num2str(i),',',num2str(idx),').*x',num2str(alfa1),'_cu','.*x',num2str(alfa2),'_cu']];
            else
                name = [name,['+2.*gxxV(',num2str(i),',',num2str(idx),').*x',num2str(alfa1),'_cu','.*x',num2str(alfa2),'_cu']];
            end
        end
    end
    text = ['    ',name,';'];
    dlmwrite(saveFile,text,'-append','delimiter','');
end
text = 'end';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'if appOrder > 2';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:length(y)
    name = [char(y(1,i)),' = ',char(y(1,i))];
    idx = 0;
    for alfa1=1:length(x)
        for alfa2=alfa1:length(x)
            for alfa3=alfa2:length(x)
                idx = idx + 1;
                if alfa1 == alfa2 && alfa2 == alfa3
                    name = [name,['+gxxxV(',num2str(i),',',num2str(idx),').*x',num2str(alfa1),'_cu','.*x',num2str(alfa2),'_cu','.*x',num2str(alfa3),'_cu']];
                elseif alfa1 == alfa2 && alfa2 ~= alfa3
                    name = [name,['+3.*gxxxV(',num2str(i),',',num2str(idx),').*x',num2str(alfa1),'_cu','.*x',num2str(alfa2),'_cu','.*x',num2str(alfa3),'_cu']];
                elseif alfa1 ~= alfa2 && alfa2 == alfa3
                    name = [name,['+3.*gxxxV(',num2str(i),',',num2str(idx),').*x',num2str(alfa1),'_cu','.*x',num2str(alfa2),'_cu','.*x',num2str(alfa3),'_cu']];
                elseif alfa1 ~= alfa2 && alfa2 ~= alfa3
                    name = [name,['+6.*gxxxV(',num2str(i),',',num2str(idx),').*x',num2str(alfa1),'_cu','.*x',num2str(alfa2),'_cu','.*x',num2str(alfa3),'_cu']];
                else
                    [alfa1,alfa2,alfa3];
                end
            end
        end
    end
    text = ['    ',name,';'];
    dlmwrite(saveFile,text,'-append','delimiter','');    
end
text = 'end';
dlmwrite(saveFile,text,'-append','delimiter','');
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'xVec_cup = zeros(numX,nx);';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'mx = setup.mx;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'if mx > 0 ';
dlmwrite(saveFile,text,'-append','delimiter','');
nx11 = length(h11);
for i=1:nx11
    TEMP = char(h11(i,:));
    TEMP = strrep(TEMP,'*','.*');
    TEMP = strrep(TEMP,'/','./');
    TEMP = strrep(TEMP,'^','.^');
    text = ['    ',char(xp(1,i)),' = ',TEMP,';'];
    dlmwrite(saveFile,text,'-append','delimiter','');
end
text = 'end ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'if setup.myx > 0 ';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:nx1-nx11
    text = ['    xVec_cup(:,mx+',num2str(i),') = ',char(y(1,i)),'(:,',num2str(i),')-gSS(',num2str(i),',1);'];
    dlmwrite(saveFile,text,'-append','delimiter','');
end
text = 'end ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'diaghx = diag(setup.hx);';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'xVec_cup(:,nx1+1:nx)  = x_cu(:,nx1+1:nx).*repmat(transpose(diaghx(nx1+1:nx)),numX,1); ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '% The controls ';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:length(y)
    name = [char(y(1,i)),' = repmat(',char(y(1,i)),',1,nNodes);'];
    text = name;
    dlmwrite(saveFile,text,'-append','delimiter','');
end
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '%Adding shocks to the states at t + 1 ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'diagEta = diag(setup.eta(nx1+1:nx,:)); ';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:nx1
    name = ['x',num2str(i),'_cup  = repmat(xVec_cup(:,',num2str(i),'),1,nNodes);'];
    text = name;
    dlmwrite(saveFile,text,'-append','delimiter','');
end
for i=1:ne
    name = ['x',num2str(nx1+i),'_cup  = repmat(xVec_cup(:,',num2str(nx1+i),'),1,nNodes) + diagEta(',num2str(i),',1).*repmat(epsNodes(',num2str(i),',:),numX,1);'];
    text = name;
    dlmwrite(saveFile,text,'-append','delimiter','');
end
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '% Scaling the states at time t+1 by the stdX';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:length(x)
    name = ['x',num2str(i),'_cup  = x',num2str(i),'_cup/scaleX_in_g(',num2str(i),',1);'];
    text = name;
    dlmwrite(saveFile,text,'-append','delimiter','');
end
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '% Getting the controls at time t+1; each element is numX * nNodes';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:length(yp)
    name = [char(yp(1,i)),' = ','g0(',num2str(i),',1)'];
    for j=1:length(x)
        name = [name,['+gx(',num2str(i),',',num2str(j),').*x',num2str(j),'_cup']];
    end
    text = [name,';'];
    dlmwrite(saveFile,text,'-append','delimiter','');
end
text = 'if appOrder > 1';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:length(yp)
    name = [char(yp(1,i)),' = ',char(yp(1,i))];
    idx = 0;
    for alfa1=1:length(x)
        for alfa2=alfa1:length(x)
            idx = idx + 1;
            if alfa1 == alfa2
                name = [name,['+gxxV(',num2str(i),',',num2str(idx),').*x',num2str(alfa1),'_cup','.*x',num2str(alfa2),'_cup']];
            else
                name = [name,['+2.*gxxV(',num2str(i),',',num2str(idx),').*x',num2str(alfa1),'_cup','.*x',num2str(alfa2),'_cup']];
            end
        end
    end
    text = ['    ',name,';'];
    dlmwrite(saveFile,text,'-append','delimiter','');
end
text = 'end';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'if appOrder > 2';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:length(y)
    name = [char(yp(1,i)),' = ',char(yp(1,i))];
    idx = 0;
    for alfa1=1:length(x)
        for alfa2=alfa1:length(x)
            for alfa3=alfa2:length(x)
                idx = idx + 1;
                if alfa1 == alfa2 && alfa2 == alfa3
                    name = [name,['+gxxxV(',num2str(i),',',num2str(idx),').*x',num2str(alfa1),'_cup','.*x',num2str(alfa2),'_cup','.*x',num2str(alfa3),'_cup']];
                elseif alfa1 == alfa2 && alfa2 ~= alfa3
                    name = [name,['+3.*gxxxV(',num2str(i),',',num2str(idx),').*x',num2str(alfa1),'_cup','.*x',num2str(alfa2),'_cup','.*x',num2str(alfa3),'_cup']];
                elseif alfa1 ~= alfa2 && alfa2 == alfa3
                    name = [name,['+3.*gxxxV(',num2str(i),',',num2str(idx),').*x',num2str(alfa1),'_cup','.*x',num2str(alfa2),'_cup','.*x',num2str(alfa3),'_cup']];
                elseif alfa1 ~= alfa2 && alfa2 ~= alfa3
                    name = [name,['+6.*gxxxV(',num2str(i),',',num2str(idx),').*x',num2str(alfa1),'_cup','.*x',num2str(alfa2),'_cup','.*x',num2str(alfa3),'_cup']];
                else
                    [alfa1,alfa2,alfa3];
                end
            end
        end
    end
    text = ['    ',name,';'];
    dlmwrite(saveFile,text,'-append','delimiter','');    
end
text = 'end';
dlmwrite(saveFile,text,'-append','delimiter','');
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '% The states at time t + 1 ';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:length(x)
     text = [char(xp(i)),' = hSS(',num2str(i),',1) + x',num2str(i),'_cup*scaleX_in_g(',num2str(i),',1);'];
     dlmwrite(saveFile,text,'-append','delimiter','');
end
text = '% The Euler errors ';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:length(y)
    TEMP = char(strucOfArrays.('f')(i,:));
    TEMP = strrep(TEMP,'*','.*');
    TEMP = strrep(TEMP,'/','./');
    TEMP = strrep(TEMP,'^','.^');
    text = ['f',num2str(i),' = ', TEMP,';'];
    dlmwrite(saveFile,text,'-append','delimiter','');    
end
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:length(y)
    text = ['f',num2str(i),' = f',num2str(i),'*wNodes;'];
    dlmwrite(saveFile,text,'-append','delimiter','');
end
text = '% Numerical integration ';
dlmwrite(saveFile,text,'-append','delimiter','');
name = ['['];
for i=1:length(y)
    if i < length(y)
        name = [name,['f',num2str(i),',']];
    else
        name = [name,['f',num2str(i)]];
    end
end
text = ['resEuler = ',name,'];'];
dlmwrite(saveFile,text,'-append','delimiter','');
%text = '% Scaling the residuals';
%dlmwrite(saveFile,text,'-append','delimiter','');
%text = 'Q   = sqrt(mean(resEuler(:).^2));';
%dlmwrite(saveFile,text,'-append','delimiter','');
text = 'res = resEuler(:);';
dlmwrite(saveFile,text,'-append','delimiter','');
%text = 'res = res/(sqrt(size(res,1))*Q^0.5);';
%dlmwrite(saveFile,text,'-append','delimiter','');
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'if nargout > 2 &&  appOrder == 1';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '%% The Jacobian if needed ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '% Current states doublicated for numerical integration';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:length(x)
    name = ['x',num2str(i),'_cu  = repmat(x_cu(:,',num2str(i),'),1,nNodes);'];
    text = name;
    dlmwrite(saveFile,text,'-append','delimiter','');
end
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
numTheta1st = size(strucOfArrays.('dF1st'),2);
for i=1:length(y)
    for j=1:numTheta1st
        if strucOfArrays.('dF1st')(i,j) ~= 0
            TEMP = char(strucOfArrays.('dF1st')(i,j));
            TEMP = strrep(TEMP,'*','.*');
            TEMP = strrep(TEMP,'/','./');
            TEMP = strrep(TEMP,'^','.^');
            for beta1=1:length(y)
                TEMP = strrep(TEMP,['g0',num2str(beta1)],['g0(',num2str(beta1),',1)']);
                TEMP = strrep(TEMP,['gSS',num2str(beta1)],['gSS(',num2str(beta1),',1)']);
                for alfa1=1:length(x)
                    TEMP = strrep(TEMP,['gx',num2str(beta1),'_',num2str(alfa1)],['gx(',num2str(beta1),',',num2str(alfa1),')']);
                end
            end
            text = ['dF1st_',num2str(i),'_',num2str(j),' = ', TEMP,';'];
            dlmwrite(saveFile,text,'-append','delimiter','');
        end
    end
end
text = '% Numerical integration ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:length(y)
    for j=1:numTheta1st
        if strucOfArrays.('dF1st')(i,j) == 0
            text = ['dF1st_',num2str(i),'_',num2str(j),' = zeros(numX,1);'];
        else
            text = ['dF1st_',num2str(i),'_',num2str(j),' = ','dF1st_',num2str(i),'_',num2str(j),'*wNodes;'];
        end
        dlmwrite(saveFile,text,'-append','delimiter','');
    end
end

for i=1:length(y)
    if i == 1
        name = '[';
    else 
        name = ' ';
    end
    for j=1:numTheta1st
        if j < numTheta1st
            name = [name,['dF1st_',num2str(i),'_',num2str(j),',']];
        else
            name = [name,['dF1st_',num2str(i),'_',num2str(j)]];
        end
    end
    if i == 1
        text = ['dF = ',name,';'];
    elseif i >1 && i < length(y)
        text = ['     ',name,';'];
    elseif i == length(y)
        text = ['     ',name,'];'];
    end
    dlmwrite(saveFile,text,'-append','delimiter','');
end
text = 'end';
dlmwrite(saveFile,text,'-append','delimiter','');


text = 'if nargout > 2 &&  appOrder == 2';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '%% The Jacobian if needed ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '% Current states doublicated for numerical integration';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:length(x)
    name = ['x',num2str(i),'_cu  = repmat(x_cu(:,',num2str(i),'),1,nNodes);'];
    text = name;
    dlmwrite(saveFile,text,'-append','delimiter','');
end
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
numTheta2nd = size(strucOfArrays.('dF2nd'),2);
for i=1:length(y)
    for j=1:numTheta2nd
        if strucOfArrays.('dF2nd')(i,j) ~= 0
            TEMP = char(strucOfArrays.('dF2nd')(i,j));
            TEMP = strrep(TEMP,'*','.*');
            TEMP = strrep(TEMP,'/','./');
            TEMP = strrep(TEMP,'^','.^');
            for beta1=1:length(y)
                TEMP = strrep(TEMP,['g0',num2str(beta1)],['g0(',num2str(beta1),',1)']);
                TEMP = strrep(TEMP,['gSS',num2str(beta1)],['gSS(',num2str(beta1),',1)']);
                idx = 0;
                for alfa1=1:length(x)
                    TEMP = strrep(TEMP,['gx',num2str(beta1),'_',num2str(alfa1)],['gx(',num2str(beta1),',',num2str(alfa1),')']);
                    for alfa2=alfa1:length(x)
                        idx = idx + 1;
                        TEMP = strrep(TEMP,['gxxV',num2str(beta1),'_',num2str(idx)],['gxxV(',num2str(beta1),',',num2str(idx),')']);
                    end
                end
            end
            text = ['dF2nd_',num2str(i),'_',num2str(j),' = ', TEMP,';'];
            dlmwrite(saveFile,text,'-append','delimiter','');
        end
    end
end
text = '% Numerical integration ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:length(y)
    for j=1:numTheta2nd
        if strucOfArrays.('dF2nd')(i,j) == 0
            text = ['dF2nd_',num2str(i),'_',num2str(j),' = zeros(numX,1);'];
        else
            text = ['dF2nd_',num2str(i),'_',num2str(j),' = ','dF2nd_',num2str(i),'_',num2str(j),'*wNodes;'];
        end
        dlmwrite(saveFile,text,'-append','delimiter','');
    end
end

for i=1:length(y)
    if i == 1
        name = '[';
    else 
        name = ' ';
    end
    for j=1:numTheta2nd
        if j < numTheta2nd
            name = [name,['dF2nd_',num2str(i),'_',num2str(j),',']];
        else
            name = [name,['dF2nd_',num2str(i),'_',num2str(j)]];
        end
    end
    if i == 1
        text = ['dF = ',name,';'];
    elseif i >1 && i < length(y)
        text = ['     ',name,';'];
    elseif i == length(y)
        text = ['     ',name,'];'];
    end
    dlmwrite(saveFile,text,'-append','delimiter','');
end
text = 'end';
dlmwrite(saveFile,text,'-append','delimiter','');




text = 'end';
dlmwrite(saveFile,text,'-append','delimiter','');
end

